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We explore the phase diagram of the dissipative Rabi-Hubbard model, as could be realized by a 
Raman-pumping scheme applied to a coupled cavity array. There exist various exotic attractors, 
including ferroelectric, antiferroelectric, and incommensurate fixed points, as well as regions of 
persistent oscillations. Many of these features can be understood analytically by truncating to the 
two lowest lying states of the Rabi model on each site. We also show that these features survive 
beyond mean-field, using Matrix Product Operator simulations. 

PACS numbers: 42.65.Sf,42.50.pq,05.70.Ln 


Introduction - A number of recent experimental break¬ 
throughs [1-4] have spurred the investigation of non¬ 
equilibrium properties of hybrid quantum many-body 
systems of interacting matter and light. Characterized 
by excitations with a finite lifetime, when sustained by 
finite-amplitude optical drives they display steady-state 
phases that are generally far richer [5-10] than their 
equilibrium counterparts [11, 12]. Critical phenomena 
in these open driven-dissipative systems often come with 
genuinely new properties and novel dynamic universal¬ 
ity classes, even when an effective temperature can be 
identified [13-17], a statement that can be made robust 
in renormalization group calculations [18, 19]. Coupled 
cavity QED systems [20-22] have emerged as natural 
platforms to study many-body physics of open quantum 
systems. The current fabrication and control capabili¬ 
ties in solid-state quantum optics allows to probe lat¬ 
tice systems [23-31] in the mesoscopic regime, providing 
a first glimpse into how macroscopic quantum behavior 
may arise far from equilibrium. It is therefore of interest 
to (i) identify a physical system where a non-equilibrium 
phase transition can be studied - at least in principle 

in the thermodynamic limit, (ii) can be compared to 
an equilibrium analogue through a proper limiting pro¬ 
cedure, and (iii) can be easily realized in an architecture 
that is currently available. 

The Rabi-Hubbard (RH) model [33] represents the 
minimal description of coupled Cavity QED systems, ex¬ 
plicitly containing terms which do not conserve particle 
number. These terms are relevant for the low-frequency 
behavior of the coupled system and their inclusion lead, 
in equilibrium, to a ^-symmetry breaking phase tran¬ 
sition between a quantum disordered para-electric phase 
and an Ising ferroelectric [33, 34]. The equilibrium RH 
transition requires a sizable inter-cavity hopping or light 
matter interaction, of the order of the transition fre¬ 
quency of cavities and qubits [33]. While such ultra¬ 



FIG. 1. Schematic array of coupled cavities in (a) ID or (b) 
2D, containing ’’Raman-driven” qubits. Inset: Cartoon of 
Raman driving: two low-lying levels of each artificial atom 
are connected via excited states. The strength of the drive 
determines the effective atom-cavity coupling [32]. 


strong coupling regimes have recently been realized in 
specific circuit QED architectures [35], they are hard to 
achieve in lattice Cavity QED settings. To overcome this 
challenge it is therefore crucial to engineer effective real¬ 
izations of the RH model by, e.g., suitably designed driv¬ 
ing schemes. In this paper we study the behavior of such 
a scheme that leads to a RH model with highly tunable 
parameters and in a fully non-equilibrium regime. 

The interplay of drive and dissipation results in ex¬ 
otic attractors, remarkably different from thermal equi¬ 
librium, with rich patterns of symmetry breaking includ¬ 
ing incommensurate and antiferroelectric ordering. In 
the following we identify and explain these orders, and 
the associated phase transitions using a variety of mean 
field approaches and then confirm the qualitative picture 
by simulating a one-dimensional open RH model with a 
Matrix Product Operator (MPO) approach [36-40]. 

Tunable Open Rabi Hubbard Model - Recently sev¬ 
eral proposals to engineer effective light-matter inter¬ 
actions by suitable designed driving schemes have ap¬ 
peared [41-44], based on a variety of platforms includ- 
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ing superconducting circuit QED, impurities in diamond, 
and ultracold-atoms [4]. Here for concreteness we con¬ 
sider a lattice of coupled cavity-QED systems, where on 
each lattice site n there is a four-level system which is 
driven and coupled to a cavity mode (see Fig. 1). We 
can write the full Hamiltonian as H = J X]{nm) + 

J2 n hn LS where J is the hopping rate, are cre¬ 

ation/annihilation operators of the cavity mode while 

h-n S = w 0 a h^n d~ X]j= 0 ,l,r,s ^'j\j)n{j\n + H lnt + H dl [ ve (t) 

describes the driven four-level atom coupled to cav¬ 
ity. The key idea of this Raman pumping scheme [41] 
is that the cavity mediates transitions between states 
0 r, 1 <->■ s (blue arrows in Fig 1), i.e. H\ nt = 
a n (<?r|f)( 0 | + g s |s)(l|) + h.c. while a two-frequency pump 
drives the transitions 1 r, 0 s (red arrows in Fig 1), 
i.e. H drive (t) = ^e-^ t |r)(l| + ^e-“? i |s)(0|+H.c.. 
The combined effect of light-matter interaction and drive 
is to induce an effective direct coupling between the two 
low lying atomic states. More formally, this can be shown 
by the standard procedure of first eliminating the explicit 
time-dependence of H dl - lve (t) moving to a rotating frame 
and then eliminating the excited states to obtain an ef¬ 
fective model for the cavity photon and the states 10 ), 11 ) 
acting as an effective qubit [32]. The effective Hamilto¬ 
nian takes the generalized RH form, 

Hrh = ~j'Yh “n®™ + X] hn W 

(nm) n 

h n = + ujaia n + (gala~ + g'ala+ + he) (2) 

where on each site n we have now a two-level system, 
<f + = |1)(0|, a~ = |0)(1|. The co-rotating and counter¬ 
rotating couplings g,g', as well as the effective cavity 
and qubit frequencies ui, u>o are tunable through the am¬ 
plitude and frequency of the Raman drive [32]. We 
stress that although described by a static Hamiltonian 
the problem retains its non-equilibrium character since 
cavity and qubit excitations are coupled to baths which 
are described, in the rotating frame, by a non-thermal 
distribution of modes. To account for the dissipative na¬ 
ture of the problem we use a master equation for the 
density matrix of the system p = — *[.Hrh, p] + J2 n ^n[p\ 
where the Liouvillian has the form 

'D n [p\ = ')C[cr-,p\ + nC[a n ,p\, (3) 

with C[X,p\ = 2XpP - ift Xp - pX tX the Lindblad 
superoperator. Here k and 7 are (constant) decay rates. 

The RH Hamiltonian in Eq. (1), as well as the dissipa- 
tor in Eq. (3), have a global Z 2 parity symmetry, corre¬ 
sponding to a simultaneous change of sign of cavity and 
qubit operators, (a’ l ',<r + ) —► (—&t, — <r + ). As a result, on 
general grounds we can expect a steady state phase dia¬ 
gram with a symmetric phase, where any quantity which 
is odd under parity will vanish, i.e. (ct“) = (a n +a' n ) = 0 , 
and a phase with broken Z 2 parity symmetry. 



FIG. 2. (a) Mean-field phase diagram of driven dissipative 
Rabi-Hubbard model a,t g = g ',uj — uiq, calculated by linear 
stability analysis of the normal state. Color scale indicates 
the wavevector of the most unstable mode. This wavevector 
predicts the ordering seen by finding the steady state solution 
for a chain of 16 cavities, as shown in panels (b-d). Panels 
(b,c) show the nearest-neighbor correlations on the vertical 
and horizontal cuts marked in panel (a). The shaded region 
shows the envelope of the limit cycle oscillations of the corre¬ 
lation function. Panel (d) shows the correlation vs separation 
at the three points marked in panel (c), revealing the incom¬ 
mensurate ordering. Parameters (also for the other figures): 
id = u >0 = 1 ,K = 0.1, 7 = 0.05. 

Mean Field Theory of Open Rabi Hubbard - To char¬ 
acterize the steady state properties of open RH model we 
make a mean-field ansatz for the system density matrix 
p = <§f) n p n . The dynamics reduces to a collection of in¬ 
equivalent single-site RH problems dtp n = —i[h n , p n \ + 
T^nlPn] + lJ\a^a n + a n a' n ,p n \ in a self-consistent field 
« n = an ansatz follows the 

standard concept of mean-field theory, that each site sees 
only the average field of its neighbors [45]. Thus, as for all 
mean-field theories, it becomes increasingly accurate in 
higher dimensions, as high coordination suppresses fluc¬ 
tuation contributions beyond the mean-field. 

We start our discussion from the g = g' case. In order 
to identify the phase boundary and to guide our anal¬ 
ysis of the ordered phase, it is useful to first study the 
instability of the homogeneous normal state, by adding 
a small perturbation to the factorized density matrix 
as done in Ref. 46, i.e. p = (£) n (p ss + Sp n ) where 
p ss is normal state density matrix obtained from the 
equation M n [p] = -i[h ni p ss ] + V n [p ss \ = 0 . Consid¬ 
ering the one-dimensional case for simplicity and tak¬ 
ing the fluctuations as plane waves of the form Sp n = 
dpke llykn ~ Vkt ^ + H.c. we obtain the equation of mo¬ 
tion: 

- w k 5pk = M n [hpk\ - t k {Tr(a<5/?fc)«[a. f ,p SJS ] + H.c.} , 

( 4 ) 
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where tk = —2 J cos(fc) is the one-dimensional bare pho¬ 
ton dispersion. A positive imaginary part of the fre¬ 
quency, Im[i/fc] > 0 , signals the growth of fluctuations 
with momentum k and the onset of normal state insta¬ 
bility. The results of linear stability analysis are plot¬ 
ted in Fig. 2(a) where we can see the phase boundary 
in the ( g , J) plane and, in color scale, the wavevector 
of the most unstable mode. Two remarkable features 
immediately appear. First, the boundary has a “nose” 
at small J, i.e. a minimal critical J required to enter 
the ordered phase. This is in contrast to the ground 
state phase diagram [33], in which the critical value of 
J asymptotically approaches 0 as g = g' — > oo. In ad¬ 
dition, the nature of the broken symmetry phase itself 
shows an interesting evolution across the phase diagram. 
As the most unstable wavevector evolves smoothly from 
k = 0, characteristic of a uniform ferroelectric (F) phase, 
toward k = 7t/2 , the wavelength must pass through ir¬ 
rational values, corresponding to an instability towards 
incommensurate order. Such symmetry-broken inhomo¬ 
geneous states requires to model a finite length array and 
we consider in panels (b)-(d) a 16 site array with periodic 
boundary conditions. We focus on correlations ,) 

which at short distance (1 = 1, panels b-c) are ferro¬ 
electric but alternate in sign as a function of distance 
l (panel d) revealing the inhomogeneous ordering. The 
finite length of the array is enough to see the trend of 
density-wave period vs hopping J, although it prevents 
a continuous evolution of the period. Such incommensu¬ 
rate order is absent in equilibrium, where the minimum 
free energy state always has a constant phase across the 
array. Another unique feature of the steady state phase 
diagram is the existence of limit cycles [47-50]. Within 
the linear stability analysis, a limit cycle can be antic¬ 
ipated if the normal state becomes unstable via a Hopf 
bifurcation [51] — i.e. if there are a pair of eigenvalues 
Vk = + iv'ff. that simultaneously become unstable, 

leading to an oscillatory instability. This in fact occurs 
for the region around g = 0.25, J > 0.8 where the most 
unstable k returns to k = 0. The existence of the limit 
cycle is confirmed by direct time evolution of the equa¬ 
tions of motion; in Fig. 2(b,c) the shaded region shows 
the envelope of the limit cycle oscillations of the correla¬ 
tion function. 

As we move away from the pure RH limit and con¬ 
sider g ^ g' two main features arise [32], namely (i) 
the shape of the phase boundary changes with multi¬ 
ple separate ordered regions and quite remarkably (ii) 
for certain values of light-matter couplings g',g we find 
an instability at k = ir corresponding to antiferroelectric 
(AF) order, where qubit and photon polarization alter¬ 
nates in sign between even and odd sites of the array, i.e. 
(&„), (a n + a'l) ~ (—l) n . This is a particularly striking 
result, considering that the effective qubit-qubit interac¬ 
tion in the equilibrium groundstate would be ferromag¬ 
netic [33], and further pinpoints the profound differences 


between the open driven-dissipative, and equilibrium in¬ 
carnations of the RH model. 


g'/g=1.0 g'/g=0.25 



FIG. 3. Properties of the effective spin model g'/g = 1 (left) 
and g'/g = 0.25 (right). (a,b) Eigenvalues and (c,d) normal 
state populations of the eigenstates of the Rabi model. The 
effective spin 1/2 model truncates to only the first two levels: 
the solid (red) and long-dashed (orange) curves. For g'/g = 
0.25, the energies and populations of these levels cross, as 
marked by an arrow. (e,f) Phase diagram of the effective Ising 
model as obtained by linear stability analysis, with color scale 
indicating the wavevector of the most unstable mode, and by 
mean field analysis (dashed line). Arrows in panel (f) mark 
the crossing points marked in panels (b,d). 

Even extending thermodynamics to negative effective 
temperature, it is not possible to explain all the features 
identified, such as limit cycles or incommensurate order. 
Predicting the pattern of steady states clearly requires 
going beyond equilibrium thermodynamics. 

Effective Spin Model - We now introduce an effec¬ 
tive spin 1/2 model which captures the essential physics 
of the RH model [32]. We start by considering the sin¬ 
gle site RH model, i.e. we set J = 0 in Eq. (1), and 
plot in figure Fig. 3 (a-d) the energies and the steady 
state populations of its eigenstates as a function of g , 
for two different values of g' /g. We first consider the 
g' = g case, panels (a,c). Here we notice that (i) the 
two low-lying states become almost degenerate for large 
g and (ii) they are the only states effectively populated. 
The idea is then to truncate the local Hilbert space to 
the two lowest energy states of the on-site RH Hamilto¬ 
nian that we denote |±) n according to their (opposite) 
parity. Within this space the on-site Hamiltonian sim¬ 
ply becomes = A t z while the Liouvillian becomes 
v f[p\ = 7 £[S-Tn + S+T+,p] + kC[A_t~ + A + T+,p\, 
where rf =x ' y ’ z are Pauli operators, and A = — E + 
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is the splitting between the lowest energy odd and even 
parity states, and A±, S± are the matrix elements A± = 
„(±|a„|=p)„, S± = „(±|(T“|=F) n - Note that the value of 
A, A±, S± are all functions of g , g ', u, ujq, found by diago¬ 
nalizing the Rabi model [32, 33]. In addition, these local 
effective two-level systems are coupled by an anisotropic 
exchange mediated by photons, J XtV ~ J(A_±A + ) which 
gives rise to an effective Hamiltonian in the transverse 
field Ising universality class. 

We now show that the effective model captures the 
salient features of the RH model. Firstly, in the limit of 
large g = g ', one can show that there is an exponentially 
small splitting A = ujq exp(—2 g 2 /uj 2 ) and the matrix el¬ 
ements become almost identical A± = (—1 ± A /u)g/ui. 
As a consequence the hopping is predominantly an Ising 
coupling T*f* +1 , and, in d dimensions, one can derive a 
simple expression for the critical hopping 


'I crit — 


«v 

W 3 


16 g 2 


( 5 ) 


Such an expression clearly explains the appearance of a 
minimum > n/2d for any finite loss, re, as opposed to 
the exponentially small critical coupling ~ A found 
in equilibrium. Furthermore Eq. (3) matches the linear- 
stability phase boundary remarkably well, see Fig. 3(e). 
In addition, it shows that as the loss re —► 0, the nose 
will move toward g — > oo, J — > 0 consistent with the 
equilibrium phase diagram [33]. 

We next consider the case g' ^ g, and plot in Fig. 3 
(b,d), the energy levels and steady state population of 
the single site RH model. In this case there are energy 
levels crossings (see arrow), corresponding to a change in 
sign of local transverse field A = E_ — E + for our ef¬ 
fective spin model. This has interesting consequences for 
the lattice model as we see in panel (f) which shows the 
phase diagram as obtained from linear stability analysis 
for g' = g/ 4. At the degeneracy point the ordered phase 
is suppressed, and beyond the crossing point an AF in¬ 
stability occurs (see bottom arrow), as recently observed 
in the transverse field Ising model [52]. Upon further in¬ 
creasing the coupling g a further transition to a normal 
phase occurs, followed by a recovery of ferroelectric or¬ 
der (top arrow). This latter effect is associated with a 
population inversion between the |±) eigenstates, as one 
can see in the level occupations shown in Fig. 3(d). 

Thus, the sequence of F-AF-F can be explained as fol¬ 
lows: At small g the ground state is that of even parity, 
and this state is the most occupied, leading to F. On in¬ 
creasing g , first the energy ordering of the odd and even 
parity states is swapped, leading to AF, where the even 
parity state is most occupied despite being of higher en¬ 
ergy. Then, the occupation of the even and odd parity 
states inverts, so that once again the lowest energy state 
is maximally occupied, and F ordering is restored. The 
qualitative picture emerging from the effective spin model 
is able to reproduce the salient features of the RH model, 



FIG. 4. Correlations of the one-dimensional effective spin 
1/2 model, evaluated in an infinite-MPO approach for g'/g = 
0.25, J = 0.5 (a) Correlations vs <?, at various separations l 
between sites, (b) Correlations vs separation at various values 
of g. These confirm the ordering seen in mean-field theory, 
specifically the sequence of F-AF-F on increasing g, but show 
only short-range order as expected in one dimension. 


both in terms of the phase boundary and in terms of pat¬ 
tern of broken symmetry phases [32], at least for mod¬ 
erate values of g. At yet higher coupling g , even higher 
states become occupied sequentially when resonances be¬ 
tween excited state energy levels occur. The occupation 
of these higher levels demonstrates the eventual failure 
of the effective spin 1/2 model at large g. 

MPO Results - A natural question is whether our 
mean field analysis survives to strong quantum fluctu¬ 
ations in low dimensions. In this respect the effective 
spin model has the advantage of being amenable to an 
exact numerical treatment with an infinite matrix prod¬ 
uct operator approach (iMPO), which we now turn to 
describe. In Fig. 4 we show iMPO results for spin corre¬ 
lators evaluated for g'/g = 0.25, J = 0.5 as a function of 
g at various separations l between sites (panel a) and as 
a function of separation l at various values of g. These 
numerically exact results confirm the ordering seen in 
mean-field, specifically the sequence of F-AF-F upon in¬ 
creasing g 1 but additionally show that in low dimensions, 
fluctuations destroy long-range order in driven dissipa¬ 
tive systems, as expected. Further results are presented 
in the supplemental material [32] supporting this state¬ 
ment. 

In summary, we have presented the steady-state phase 
diagram of the non-equilibrium Rabi-Hubbard model, 
using various mean-field-based techniques and a ma¬ 
trix product operator approach that can capture effects 
beyond mean-field. The phase diagram of the non¬ 
equilibrium model was found to be far richer than the 
equilibrium analogue, exhibiting ferroelectric, antiferro- 
electric and incommensurate ordering. In addition, the 
phase diagram was found to also exhibit limit-cycle solu¬ 
tions. The MPO results confirm qualitatively the pattern 
of the phases. 

The research data supporting this publication can be 
accessed at http://dx.doi.org/10.17630/a8829c09-7d49- 
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The supplementary material is organized as follows: 
Section I presents further details of the possible ex¬ 
plicit driving scheme realizing a (tunable) Rabi-Hubbard 
model as presented in Eq. (1) of the main text. Section II 
presents the mean field and linear stability results for the 
generalized Rabi-Hubbard, with g' ^ g that were men¬ 
tioned in the main text. Section III discusses the con¬ 
struction of the effective spin model introduced in the 
main text and the linear stability analysis of this model. 
Section IV presents further results on the effective spin 
model obtained by MPO simulations of an infinite chain. 


I. RAMAN DRIVING SCHEME AND TUNABLE 
RABI HUBBARD 

In this section we present a possible scheme to obtain 
the Rabi Hubbard model starting from a driven four- 
level atom in a cavity. This scheme follows closely the 
original proposal [1] for a driven realisation of the gen¬ 
eralized Dicke model, which was recently experimentally 
realized [2]. To show most clearly how Rabi-like interac¬ 
tions can be generated from multi-atom driven problem 
we focus here on the isolated single cavity Hamiltonian, 
and disregard for now the photon hopping to other neigh¬ 
boring resonators. The interference between coupling to 
other cavities and pump induced terms can generate very 
weak long-ranged hopping and interactions between cav¬ 
ities. The analysis of these small corrections is left for 
future studies. 

Our starting point is thus a four-level (artificial) atom 
in an optical cavity (or microwave resonator), supporting 
a single photon mode, see Figure 1. The full local Hamil¬ 
tonian is made of several terms. The non-interacting 
atom and cavity are described by 

Ho=w cav d t a+ "^2 E o\j){j\- 

j= 0,1,r,s 

Following Ref. 1 we consider a case where the cavity 
mode drives transitions between the states |0) -H- | r) and 
11) •<->• | s ). Selecting these transitions can either be done 
by tuning the other transitions to be far off resonance 
with the cavity mode, or ideally, by engineering a sys¬ 
tem with selection rules restricting which transitions the 
cavity mode polarisation couples to [2], The resulting 


cavity-atom coupling Hamiltonian reads: 

Hint = 9r | r)(0| a + g s \ s)( 1| a + H.c. 


w 

Es 7777 A s | r ) 



FIG. 1. Cartoon of the level structure and Raman driving 
scheme to engineer a tunable Rabi Hamiltonian. 

Since the goal is to generate effective light-matter in¬ 
teractions between the photon and the low-lying doublet 
| 0), 11), we need to connect these states through resonant 
two-photon transitions. This can be done by adding the 
drive term: 

Hdrive (i) = ^ ‘ I r) (11 + 3 s e-^ ‘ | s) (0| + H.c. 

A cartoon of the level scheme, and the pump- and cavity- 
mediated transitions is shown in Figure 1. 

The crucial feature of this scheme is that the laser 
drives transitions between |1) -H- | r) and |0) O | s) so 
that the combined effect of H; nt and H^me is to induce 
two-photon transitions between | 0) and 11), involving 
emission or absorption of a cavity photon. These will 
become the familiar rotating and counter-rotating terms 
of the Rabi model. 
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We proceed in two steps, that we briefly outline here. 
First, we eliminate the explicit time dependence of the 
drive by performing an unitary transformation to a ro¬ 
tating (comoving) frame. This can be implemented by 
the operator: 


Cl(t) = exp , K = a a) a + /3j| i)(i\, 

i 


drive as follows: 


LO = UJ C 


_3£ 

V 

i /n 


9t_ 

At 


u: 


- U!' 


<*>=*-* 




■y y _3r_ Us 

~ 2 V At At 


9 = 


[h-l r 

2A. r 


9 = 


9s^s 

2A7 


which has the property that, fT(i) aCl(t) = ae~ lat and 
flt(t) | i)( j| fi(t) = \i)(j\e l (P i ~PA t . We can fix the 
parameters in K in such a way that the transformed 
Hamiltonian H = idt ^fV(f)^ £l(t)+£tf (t) H(t) f)(i), with 

H = H 0 + H ln t + H d r ive, becomes time-independent. We 
stress that while the problem in this new rotating frame 
is now fully time-independent, it remains intrinsically out 
of equilibrium in nature since (i) we are forced to study 
the dynamics of H , and (ii) the distribution function of 
any bath modes the system couples to becomes strongly 
non-thermal (breaking detailed balance between gain and 
loss), as the unitary transform II must also be applied to 
the system-bath coupling. 

We now proceed to the elimination of atomic excited 
states | r ), | s) in order to obtain an Hamiltonian for the 
manifold | i) = | 0), 11) that will play the role of our qubit. 
This can be done perturbatively, using a Sclrrieffer-Wolf 
transformation, Hg = e lS He~ lS with a generator S cou¬ 
pling the qubit manifold |0),|1) to the excited states 
\r),\s) 

S =i£t/w + k' sm+ 

+ |^oH(0| + ||»l»)( 1 l) +H.c.. 

This is obtained by imposing the condition [Hq,iS] = 
Hint + -ffdrive (in terms of the time-independent Hamilto¬ 
nian). We have introduced the detunings 

An = W C av + Ef) — E r , 

A“ = w cav + Ei — E a + {ojp — Wp)/2, 

A r = Ei — E r + u>p, A s = E 0 — E s + <jJp- 

To leading order in the strength of the drive 
and light-matter coupling we get Hg = i?R a bi + 
0(,9r /A r , g s j A s , /A r , A s ) where 

^Rabi = w a) a + -^"d z + A a)aa z + g{a)a~ + aa + ) 

+ g’ (a 1 ' <x + + aa~). (1) 

This is a generalized Rabi model, for which the parame¬ 
ters depend on the frequency and strength of the external 


where we defined 1/A r / S = 1/A“j s + 1/A r / S . 

To complete the mapping to the generalized Rabi 
model we must choose parameters such that the coef¬ 
ficient A = 0, as also discussed in [1], To this end we 
note that AJf is drive-independent and therefore fixed by 
the physical realization, while A" can be tuned by uj s p . 
We therefore impose the condition: 

q 2 

AC O A UJ JS 

AAg 9 5 

9r 

which immediately gives A = 0. From this we see that we 
may want to have g s /g r ~ 1 in order to have both A“, A“ 
large. The second drive frequency uj r p can then be used to 
control the effective detuning between the qubit and the 
resonator. In addition to changing the detuning through 
the pump frequencies, we can vary the pump strengths 
to tune (independently) the relative strengths of 
the rotating and counter-rotating terms. 


II. OPEN RABI HUBBARD FOR g' =£ g 

In this section we present results, anticipated in the 
main text, for the phase diagram of the generalized open 
Rabi Hubbard, with g' ^ g. We use again linear stabil¬ 
ity analysis of the normal phase, as described in the main 
text, and plot in figure 2 the phase boundary in the ( g, J) 
plane and, in color scale, the wavevector of the most 
unstable mode, for two different values of g'/g, respec¬ 
tively g' = 0.5 g (left panel) and g' = 0.25 g (right panel). 
These results show several interesting features. Firstly 
there is a change in the topography of the phase bound¬ 
ary, with multiple separate ordered regions, as opposed 
to the case g' = g which features a single, connected, 
broken symmetry phase. As a result, upon increasing 
the light-matter interaction g one can have a sequence 
of transitions, where the symmetry is broken first, then 
restored and then broken again. In addition, for smaller 
values of the ratio g'/g the nature of the ordered phase 
changes qualitatively and an instability at k = n, even¬ 
tually emerges. This corresponds to an antiferroelectric 
(AF) order, where photon and two-level system (qubit) 
polarization changes sign between neighboring sites. This 
latter result is particularly striking if interpreted in terms 
of equilibrium physics. Indeed the effective qubit-qubit 
interaction mediated by a population of photon modes 
in equilibrium at zero temperature would be negative at 
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FIG. 2. Mean field phase diagram of the generalized (</ A g) 
open Rabi Hubbard model as obtained from the linear stabil¬ 
ity analysis of the normal phase. We plot the phase boundary 
in the (g, J ) plane and, in color scale, the wavevector of the 
most unstable mode, for two different values of g'/g , respec¬ 
tively g' = 0.5gr (left panel) and g' = 0.25 g (right panel). 


short distance [3] leading to a uniform ferroelectric pat¬ 
tern in the ground state. The results of figure 2 show the 
profound differences between the open driven-dissipative 
and equilibrium incarnations of the Rabi Hubbard model. 
As we have discussed in the main text, these intriguing 
results can be understood qualitatively in terms of an 
effective open-dissipative spin model whose construction 
we are now going to discuss in detail in the rest of this 
supplementary material. 


III. EFFECTIVE SPIN MODEL 

By truncating the on-site Rabi model (L/rubi obtained 
in the previous section) to the lowest two eigenstates and 
labelling these as |±)„ (with corresponding energies E±) 
to denote their different parities, we obtain the following 
low-frequency effective Master equation describing the 
dynamics of the generalized Rabi-Hubbard model: 

d t p = -i[H eS , p] + k£[A_t~ + A + r+,p} + 

+ ~/C[S-f- + S+T+,p\, (2) 

where the effective Hamiltonian for the generalized Rabi- 
Hubbard Model reads (for a detailed discussion of this 
projection, see Refs. 3 and 4 which discuss the case g = 
9 ') 

H ‘ 39 = - y ff 
2 ^ " 
n 

- J ( A ~ X n + A +Tn ) { A -Tm + A + x +) , 

(n,m) 

where r' n are Pauli matrices on the nth site, A = E + —E_ 
is the lowest excitation energy of the local Rabi model, 
and we have introduced the real-valued coefficients 

A± =„(±|a„|=F) n , S± =„(±|o-“|T)n- (3) 

For later convenience it is useful to introduce the even 
and odd combinations of these coefficients, namely 

a± = A + ± A-, <;± = S + ± S-. (4) 


The effective-model parameters A, a±, <r± are functions of 
the parameters g,g\ ojq in the Rabi model. We will dis¬ 
cuss below analytic approximations in the limit of large 
g = g ', the effective model is however useful more gener¬ 
ally. 


A. Mean-field theory, linear stability 


The effective model can be used to find a closed-form 
expression for the critical hopping required to stabilize 
the ordered state within the mean-field theory, as dis¬ 
cussed next. Within a mean-field decoupling of the hop¬ 
ping, the single site Hamiltonian has the form 


h: 


eflt.MF _ 


2 ^~ n 4 


--a 2 _t y h y n , (5) 


where the local effective fields 


h n — ^ 1 Xmi hn — y ' Urn: > 

m:{mn) m:(mn) 

are written in terms of x n = (f^),y n = {r y ),z n = (t*), 
account for the effect of the sites m connected by photon 
hopping to sites n. By introducing the damping rates 

r± = na ± -I- 7 ? ± , 


and the normal state inversion 

z o = r+ (ra+a- + 7^-0 , 

we can write down the mean field dynamics in a compact 
form as 

^Un r_X n — Z n , 

^tVn — r -)-Vn h n Zm 

d t z n = (r+ + r_) (z 0 - z n ) - J -olL hn y n + ^ a 2 _ h v n x n . 

The normal state solution is given by z n = z 0 and 
x n = Un = 0. Considering small fluctuations about this 
stationary solution we obtain a decoupled equation for 
z n , with damping rate T + + T_ while for the other com¬ 
ponents the ansatz x n = J2k x k e^ k " n ~ l ' kt > gives the sec¬ 
ular equation for the frequencies 

(iv k — T_){iv k — T + ) + (A - t k a 2 + z 0 ) (A - t k a 2 _z 0 ) = 0, 

where t k = — 2 J cos(fcj) is the d-dimensional bare 

photon dispersion. The instability, corresponding to a 
pitchfork bifurcation, is given by v k = 0. This leads to 
a simple expression for the critical J in the limit |a_| <C 
|a + |, which, as we will see, is satisfied at large g = g'. 
For the case where a ferromagnetic (anti-ferromagnetic) 
instability ki = 0(7r) occurs, this expression is: 


Jcnt ~ T 2da 2 + z 0 



A 


( 6 ) 
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Since we consider only positive J , it is clear from the form 
of this equation that the sign of tk required, i.e. whether 
the instability is ferro (F) or antiferro (AF), is determined 
by the sign of the product Azq . For g = g'. we always 
have A > 0,2 O < 0 and so only the ferromagnetic case 
(negative sign) occurs Thus, as discussed in the main 
text, the level crossing at A = 0, and the inversion point 
at Zo = 0 lead to suppression of ordering, and a switch 
between F and AF ordering. 


B. Large g limit and Effective Model parameters 


In the limit of large g = g' one may use an approximate 
solution of the on-site Rabi model to derive simple ana¬ 
lytic expressions for the lowest excitation energy A and 
the matrix elements a± ,g±. To extract these expressions 
we start from the single site Rabi Hamiltonian 

Rabi w 0 


H 


= + uia* a + g{a* + a)a 


(7) 


and perform the unitary transformation 

9 '* t 


U = exp 


-Xa* 


X = — (cr — a), 
uj 


to obtain FT Rabl = H Kahl U in the form 

_ 2 


Rabi _ _£_ 


9_ 

OJ 


■ uia) a - 


Wo 

2 L 


cosh 2X a z — i sinh 2X a y 


(8) 

In absence of the term in brackets the spectrum has an 
infinite sequence of two-fold degenerate states with ener¬ 
gies K = ~(g 2 / u>) + nui and corresponding eigenstates 
|er = {f, 4-},n) in the transformed basis. In what fol¬ 
lows we carry out a perturbation expansion in SH = 


Wo 

2 


cosh 2 A' a z — i sinh 2X a v 


for the states in the low¬ 


est manifold (identified by n = 0). This is justified by the 
smallness of the parameter (0|e ±2 ' Y |0) = exp (—2(g/w) 2 ) 
in the large-g limit j/w> 1. To lowest order, this splits 
the n = 0 doublet and yields an analytic expression for 
the lowest excitation energy of the Rabi model: 


A = u) 0 exp (—2 g 2 /u> 2 ) . 


(9) 


To obtain an analytic expression for the matrix elements 
Eq. 4, we also need to compute the first order correction 
to the wavefunctions | cr, 0): 


|o\0) = |cr, 0) - 


E 

n>0,T 


|t, n)(r, n\6H\a, 0) 


( 10 ) 


where |cr, 0) denotes the corrected wavefunction, and the 
matrix elements of the perturbation Hamiltonian can be 
found using the following expression: 

(n|<j.Hj0) = ~^{ n \ cosh(2A)<r z — isinh(2X)(7y|0) 


A (2g 

2\fn\ V w 


[In 


— i I, 


n(zodd Gyl 


(ii) 


where I n ^A is an indicator function, taking value 1 for 
n € A and zero elsewhere. 

We can now evaluate the matrix elements of photon 
and spin operators in the space spanned by the low- 
energy doublet, that in the transformed basis are simply 
the states | t>0), | j.,0). To evaluate the spin matrix el¬ 
ements we note that, under the unitary transformation 
we have 


U*&~U= - |d x + cosh(2X)( 
from which we conclude that 

s+ = <t,0|tn*-fr|4,0> = 
S- = {i,0\U^a~U\ f, 0) = 


iby) + sinh(2A')cr z , 




(12) 


(13) 

Wo / 


A \ 



(14) 


where to leading order we have only accounted for the 
transformation of the operator. Turning to the pho¬ 
ton matrix elements we first note that under the unitary 
transformation we have 

U^aU = a — —b x , (15) 


therefore to evaluate the matrix element one has to use 
the perturbed eigenstate to obtain, at leading order, 


A + = (t,0|a— -<7 X |4.,0) 

UJ 


gA g _ gA g 

2 ’ ’ -/i— 2 

UJ Z UJ UJ Z LJ 


(16) 

Using these results, we find that the matrix elements 
can be found to leading order to be a+ = 2 g/co,a^ = 
2gA/w 2 ,? + = l,c_ = A/ujq, from which we can ob¬ 
tain the approximate expression for the critical boundary 
given in the main text. 


IV. MPO RESULTS FOR THE RABI HUBBARD 

CASE, g = g r 

In order to test the predictions of mean-field theory, 
we may use a infinite Matrix-Product Operator (iMPO) 
approach [5-8] to find the non-equilibrium steady state of 
the effective model, supplementary Eq. (2). The iMPO 
approach is applicable to a translationally invariant prob¬ 
lem [9], allowing one to evolve only a representative pair 
of sites, and has no finite-size effects. Despite this, as 
seen in Fig. 5 of the paper, iMPO can fully describe inho¬ 
mogeneous order; this is because multi-site expectations 
involve traces of products of the representative matrix. 
Our approach, as in previous work [10], is to simulate 
an adiabatic sweep, starting from the point g = g’ = 0 
where the state is trivially a product state. Such a sweep 
avoids the high transient entanglement that occurs with 
a quench of parameters, and means that the correlation 
functions we measure are converged for the bond dimen¬ 
sion x = 120 we use. 
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FIG. 3. Correlations of the one-dimensional effective spin 
1/2 model, evaluated in an infinte-MPO approach for g' = 
g wt 1, J = 0.5 (a) Correlations vs g, at various separations l 
between sites, (b) Correlations vs separation at various values 
of g. 


In the main text we showed results for g'/g = 0.25, il¬ 
lustrating that the F-AF-F sequence predicted by mean- 
field theory still exists (albeit short-ranged) for the MPO 
simulation. Here we present additional MPO results on 
the effective spin model in the pure Rabi case, i.e. for 
g = g' = 1. In particular the data reported in fig¬ 
ure 3 show short-ranged ferroelectric correlations devel¬ 
oping for several values of the light-matter interaction, 
consistently with the mean field picture. 
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